#include "fintrf.h"
C
#if 0
C     
C     
#endif
	SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
	IMPLICIT NONE
	MWPOINTER PLHS(*), PRHS(*)
	INTEGER MXGETM, MXGETN
	INTEGER NLHS, NRHS
	MWPOINTER ROWS, COLS
	MWSIZE SIZE
	
C       Check for proper number of arguments. 
	IF (NRHS .NE. 1) THEN
	   CALL MEXERRMSGTXT('One input required in wofz.')
	ENDIF
	ROWS = MXGETM(PRHS(1))
	COLS = MXGETN(PRHS(1))
	SIZE = ROWS*COLS
!	CALL INTERMEDIATE GATEWAY:
	CALL GATEWAY(PLHS,PRHS,SIZE)
	
	END 
!	==============================================================
	
	SUBROUTINE GATEWAY(PLHS, PRHS, SIZE)
	IMPLICIT NONE
	MWPOINTER PLHS, PRHS
	MWSIZE SIZE
	MWPOINTER MXCREATEDOUBLEMATRIX, MXGETPR, MXGETPI
	MWPOINTER MXGETM, MXGETN, MXISCOMPLEX
	MWPOINTER ROWS, COLS
	INTEGER I, ISCOMPLEX
	DOUBLE PRECISION A(SIZE), B(SIZE), C(SIZE), D(SIZE)
	LOGICAL FLAG   
	
	CALL MXCOPYPTRTOREAL8(MXGETPR(PRHS),A,SIZE)
!       Checks if there are imaginary numbers in data
	ISCOMPLEX = MXISCOMPLEX(PRHS)               
	IF (ISCOMPLEX.EQ.1) THEN   
	   CALL MXCOPYPTRTOREAL8(MXGETPI(PRHS),B,SIZE)
	END IF               
	ROWS = MXGETM(PRHS)
	COLS = MXGETN(PRHS)
	PLHS = MXCREATEDOUBLEMATRIX(ROWS, COLS, 1)

!       CALL COMPUTATIONAL ROUTINE:
	DO 20, I=1,SIZE
	   IF (ISCOMPLEX.EQ.0) THEN   
	      B(I) = 0
	   END IF   
	   CALL WOFZ(A(I),B(I),C(I),D(I),FLAG)
 20	CONTINUE
	CALL MXCOPYREAL8TOPTR(C,MXGETPR(PLHS),SIZE)
	CALL MXCOPYREAL8TOPTR(D,MXGETPI(PLHS),SIZE)
	END
!	===============================================================
	SUBROUTINE WOFZ (X, Y, U, V)
	IMPLICIT DOUBLE PRECISION (A-H, O-Z)
	DOUBLE PRECISION A0, A1, A2, A3, A4, A5, A6
	DOUBLE PRECISION B0, B1, B2, B3, B4, B5, B6
	COMPLEX *16  T, W
	A0 = 122.607931777104326
	A1 = 214.382388694706425 
	A2 = 181.928533092181549 
	A3 = 93.155580458138441 
	A4 = 30.180142196210589 
	A5 = 5.912626209773153 
	A6 = 0.564189583562615 
	B0 = 122.607931773875350
	B1 = 352.730625110963558
	B2 = 457.334478783897737
	B3 = 348.703917719495792
	B4 = 170.354001821091472
	B5 = 53.992906912940207
	B6 = 10.479857114260399
	T  = DCMPLX(Y,-X)
	W  = ((((((A6*T + A5)*T + A4)*T + A3)*T + A2)*T + A1)*T+A0)
	1      /(((((((T + B6)*T + B5)*T + B4)*T + B3)*T + B2)*T + B1)*T+B0)
	U = DREAL(W)
	V = DIMAG(W)
	RETURN
	END
